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Abstract 

A very efficient large-order perturbation theory is formulated for the nuclear 
motion of a linear triatomic molecule. To demonstrate the method, all of the 
experimentally observed rotational energies, with values of J almost up to 
100, for the ground and first excited vibrational states of CO2 and for the 
ground vibrational states of N2O and of OCS are calculated. All coupling 
between vibration and rotation is included. The perturbation expansions re- 
ported here are rapidly convergent. The perturbation parameter is D~ 1 / 2 , 
where D is the dimensionality of space. Increasing D is qualitatively similar 
to increasing the angular momentum quantum number J. Therefore, this 

approach is especially suited for states with high rotational excitation. The 

5/3 

computational cost of the method scales only as JN V , where N v is the size 
of the vibrational basis set. 



I. INTRODUCTION 

The calculation of highly excited rovibrational states of molecules is a challenging com- 
putational problem. The traditional method is numerical diagonalization of the Hamiltonian 
in a finite basis set. This works well for low- lying states, but for highly excited states it 
is tractable only with judicious choice of basis set and the use of special numerical 
techniques |3],|4]||, on account of the extremely large matrices that are needed in order to 
obtain sufficient accuracy. 

It has been suggested that large-order perturbation theory is a viable alternative ap- 
proach to calculating vibrational spectra. We recently compared numerical diagonalization 
with large-order perturbation theory in powers of the coupling constant, for a model system 
of coupled anharmonic oscillators, and found that the perturbation theory had a significantly 
lower computational cost [0]. Here we apply dimensional perturbation theory (DPT) to the 
calculation of vibration-rotation spectra for actual triatomic molecules. This is a type of 
perturbation theory that is especially suited for states with high angular momentum ||. 

This theory [p|, p]0| , p] , p^ ,^| is one example of a class of semiclassical perturbation expan- 
sions that have been applied to problems in high-energy physics, statistical mechanics, solid 
state physics, and polymer science. There has been much recent interest in using DPT to 
solve the electronic Schrodinger equation for atoms. At the zeroth order of the perturbation 



theory an atom becomes a rigid rotor, with the electrons stationary relative to each other. 
At first order the electrons vibrate harmonically. The first-order results for electronic ener- 
gies of atoms are typically accurate to within 1% jjTJJ, and computational techniques have 



now been developed that make it possible to carry the perturbation theory to much higher 
order |5^|7|. 

Thus far there have been few attempts to apply this theory to molecules, despite the 
resemblance between DPT and the standard qualitative picture of molecular vibrations. 



First-order expansions have been calculated for the electronic energy of H 2 [fl8[], for the 
full electronic-vibrational problem (without the Born-Oppenheimer approximation) for n 2 + 
T9| , and for the vibrational motion of atomic clusters ||20|| . These studies have provided 



intriguing qualitative models of correlated many-particle dynamics and semiquantitative 
results for energy eigenvalues. Recently, however, Kais and Herschbach [H] used DPT to 
compute the rotation spectrum of the ground vibrational state of H 2 . The accuracy of their 
results, from just first-order perturbation theory, remained excellent even for J as high as 
50. 

The distinctive feature of DPT is that it uses a Hamiltonian operator that has been 
analytically continued to a nonphysical hyperdimensional coordinate space. Typically, the 
kinetic energy operator is generalized to a D-dimensional form while the potential energy is 
left unchanged. Since this dimensional continuation is nonphysical, it allows for a great deal 
of flexibility in how one chooses to define it. We will formulate the perturbation theory in 
such a way that it closely resembles the Kais- Herschbach diatomic analysis. In Section [n] we 
develop our dimensional continuation of the kinetic energy operator. Section |T| describes 
the method we have used to compute the expansion coefficients. In Section [TV] we use the 
theory to reproduce the empirical rotational spectra for various linear triatomic molecules 
and in Section [V] we discuss advantages of this approach. 



II. DIMENSIONAL CONTINUATION OF THE KINETIC ENERGY OPERATOR 

DPT is based on the fact that the dynamics of a system of particles simplifies as the 
dimensionality becomes very large. This simplication depends on the partitioning of the 
coordinates into translational, rotational, and internal degrees of freedom. Consider an 
arbitrary system of iV particles in D spatial dimensions and let r trans , r rot , and r int be the 
number of coordinates of the indicated type. r trans and r rot increase linearly with D while r int 
quickly reaches a maximum value of N(N — l)/2 and then remains constant with increasing 
D. Since the potential energy depends only on the internal coordinates, the energy levels of 
the system can in principle be determined from a differential equation in just r- mt coordinates 
even in the limit D — > oo. 

The separation of the D translational coordinates is straightforward. Consider a system 
of 3 atoms, which we will describe in terms of the center-of-mass vector, r cm ; the vector 
between atoms 1 and 2, r s ; and the vector r a that leads from the center of mass of atoms 1 
and 2 to particle 3. If ri, r 2 and r3 specify the positions of the atoms relative to the center 
of mass, the r s = r 2 — i"i and 

r a = r 3 - [(m 1 /m s )r 1 + (m 2 /m s )r 2 ], (2.1) 

where m s = mi + m 2 . The the kinetic energy operator can be written as 
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T = -iA Rcm + T rel , T, 
in terms of the mass-normalized coordinates 
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where m = m\ + m 2 + 777,3 is the total mass of the system. 

The separation of rotational coordinates is more complicated, but is key to the simplying 
effect of increasing D. The qualitative effect of separating out angular coordinates is to 
introduce a centrifugal potential into the internal-coordinate kinetic energy operator. The 
linear increase with D in the number of angular degrees of freedom causes the strength of 
the centrifugal potential to increase quadratically with D. The derivatives with respect to 
the internal coordinates are overwhelmed by the centrifugal potential in the limit D — > oo. 
This elimination of the derivatives yields a Schrodinger equation that can be solved exactly 
even for many-particle systems. 

Consider first a diatomic molecule with mass-normalized internuclear distance R s . If we 
define the D-dimensional Laplacian as 



d 2 

^ dx 2 ' 

i=l,D ux i. 



(2.4) 



expressed as [21,22 



where the Xi are the Cartesian components of R s , then the Schrodinger equation can be 

$(R*)Y JM , D -i(to) = 0, 



.-A Rs + V(R s )-E 



where 



1 d 
RP- 1 dIL 
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dR, 



R 2 



(2.5) 



(2.6) 



The Yj t M,D-i, which are hyperdimensional generalizations of the spherical harmonics, are 
functions of the set O of D — 1 independent angular coordinates. The substitution 



$(i? s ) = /C (D ~ 1)/2 *(#s 
removes the first derivatives, with the result 

1 d 2 , J D (J D + 1) 



2dRl 



+ 



2R 2 



V -E 



* = 0, 



(2.7) 



(2.8) 



where 



J D = J+(D-3)/2. 



(2.9) 
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This is the standard approach to dimensional continuation |23 , |24| in which the kinetic energy 
is defined in a D-dimensional Cartesian coordinate system while the potential energy is 
left in its 3-dimensional form (except, perhaps, for certain scale factors p4]j25| , |26[| ) . It is 
straightforward to use Eq. ( |2.8|) to calculate a large-order asymptotic expansion for E in 
powers of l/D or l/J D |1J,[T§. 

For triatomics this analysis is complicated by the fact that the total angular momentum 
of the system depends not just on the external rotation coordinates but also on the angle be- 
tween R s and R a . Rather than carry out the separation of variables with the .D-dimensional 
Laplacian, we will use a simpler approach. There are only two considerations that restrict 
the way in which the dimensional continuation is defined: The Schrodinger equation must 
be the correct physical equation for D = 3, and the large-D limit must be a reasonable, 
and solvable, zeroth-order approximation for the perturbation theory. Our strategy is to 
perform the separation of variables using the 3-dimensional Laplacian and then introduce 
an arbitrary D dependence consistent with these considerations. In essence, we will allow 
one body-fixed axis to rotate in D dimensions while a second body-fixed axis rotates about 
the first in a 3-dimensional space. 

For the internal coordinates we will use i? s and the cylindrical coordinates p and z, 



2 i 2 
P + Z 



Rl, p/z = tan (5 



(2.10) 



where f3 is the angle between R a and R s . The derivation of the 3-dimensional kinetic energy 
operator in these coordinates is given in Appendix A. Let L be the total angular momentum 
operator, L\ the operator for the projection of angular momentum onto a space-fixed axis 
xi, and K the operator for the projection onto the body-fixed axis R s , and let J, M, and 
K be the corresponding quantum numbers. It is convenient to describe the wavefunction in 
terms of symmetrized basis functions 



c(±) 



-1/2 



>J,M.K 



± 



>J,M,-K) 



(2.11) 



where the 4>j,m,k(R s , P, z) are simultaneous eigenfunctions of L, L Xl , and L s , and K is now 
defined to be nonnegative. Let a be a parity index that indicates the sign in Eq. ( [2.11 ). 
Then the internal-coordinate part of the kinetic energy operator for given J and M can be 
expressed in terms of the matrix elements Tj^k^k 1 ,^ — (K, a\T re \\K', a'). In Appendix A we 
derive the following expression: 
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(4,±) 
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TJJ(J+1)-K(K±1] 



X 



(2.13c) 

(2.13d) 
(2.13e) 
(2.13f) 

(2.13g) 



with Eq. ( |2.13g| ) valid for K > 1 and for the case K — 1, a — +. Note that there will be two 
independent sets of (pjMKi on account of the form of the Kronecker deltas that multiply the 

TjkI in Eq- (|2.12|) . One set will consist of basis functions for which a = + and K is even 
or a = — and K is odd, while the other will consist of those functions with a = — and K 
even or a = + and K odd. The lowering operator Tj\~J is zero for a = — but nonzero for 

a = + because of the fact that ^jmo = 0- 

The only modification we will make in the 3-dimensional operator TjK, a is to replace 
T^jk with the dimensional continuation 
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which is analogous to the diatomic result in Eq. (2J3). If we use the scaled basis functions 

^tlu(^P^) = Ri D ~ 1)/2 M, K (Rs,P,z), (2.15) 



then Tj 1 ^ becomes 
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with J D given by Eq. (|2.9| ). 



III. THE PERTURBATION THEORY 



A. Dimensional Scaling 

We consider here systems that can be described by a potential of the form 



tt 2 / \2 i 1 2 2 

+ l - u 2 s (R s - R cq ) 2 +X(Rs- R eq ) p\ (3.1) 

which is appropriate for a linear triatomic molecule. We will assume that U is independent 
of dimension. 

Our goal is to derive an asymptotic expansion about the limit D — > oo in terms of the 
parameter 



K 



D- 1/2 , (3.2) 



which will be summed at the physical value k = 3 1//2 . To obtain a useful k — > limit, we 
introduce the dimension-scaled quantities 



r, 



/ci2eq, A = k^X, E = k 2 E, (3.3) 



and the displacement coordinates 

r = K(R s -R eq ), y = z-z eq . (3.4) 

We will treat r eq , A and E as ^-independent constants with values 3~ 1//2 i? e q, 3 1//2 A, and 3 _1 i?, 
respectively. 

The dimension-scaled Hamiltonian matrix can be expressed in terms of the operator 



7~C — &o'a$K'K 



1 

2(r cq + r) 

+ l^^a',—a (^K',K+lTj i K,l+ , (3.5) 



where 7i( diag ) is the separable operator 



o/(diag) _ K 2 T (2) 4 1 d 2 

1 1 

+ _ UJ y + K 2 -(uly 2 + L;y) 

1 + 4( J - l)*: 2 + 4(J 2 - K 2 -2J- l/4)^ 4 
8(r eq + r) 2 

(3.6) 

In the limit k — > all derivatives in Eq. ( |3.5|) drop out, leaving the system localized 
at p — 0, y — 0, and r = r min , where r min is the minimum of the effective potential 
(l/8)(r eq + r)~ 2 + (l/2)u; 2 r 2 . The numerical value of r min can be determined from the equa- 
tion 

4w s 2 r min (r cq + r min ) 3 -1 = 0. (3.7) 
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B. The Energy Expansion 



The next step is to introduce a dimension-scaled displacement coordinate x, 

r = r min + kx (3.8) 



and expand Eq. (|3.5|) in powers of k, 



H = E + k 2 J2 «*ft fc , (3-9) 

fc=0 

where E is the zeroth-order solution for E, 

E = \B + §u; s r^ in , (3.10) 
in terms of the effective rotational constant 

B = — -. (3.11) 

O r 4- r ■ ) 2 

eq ~ 1 mm J 

At order n 2 we have a separable Schrodinger equation that can be solved exactly. Our 
strategy will be to express the wavefunction at any order in k in the form of a linear 
combination of the separable product eigenfunctions obtained at order k 2 . 

Note that several of the terms in Eqs. (|3.5|) and (|3.6|) are proportional to k 4 . We can im- 
prove the convergence of our perturbation theory by redefining the dimensional continuation 
so that these terms enter at lower orders [p7[] . Consider the operator 7^ < - diag - ) . The second 



derivative with respect to r, after the change of variable in Eq. ( |3.8| ), becomes proportional 
to k, 2 . We will replace the other factor of k 4 in Eq. ( |3.6| ) with n 2 /3, taking advantage of 
the fact that our Hamiltonian will still be valid as long as it is correct for the physical case 
k = 3 -1 / 2 . In Eq. ( |3.5| ) we will replace the factors of n A with 3 -1 / 2 /? 3 , in order to still have 
a separable solution at order k 2 . 

Thus, at order k 2 , we have the operator 



where 



H = 5 a ,,J K , K (% + Vo), (3.12) 



_!_(&_ & 2 _\_l_(l^^_K 2 \ 
10 ~ 2 Ux 2 + dy* ) 2 \pdp P dp p 2 ) {6 - L6) 



and 



V = +^u 2 x 2 + X -uly 2 + ij 2 p 2 

+i[J(J+l)-(K 2 + 13/4)] J B, (3.14) 
with the effective frequencies u s and (I>b given by 

u 2 = lu 2 s +3B 2 , Col = + 2Ar min . (3.15) 
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The normalized eigenfunctions of 7io are 

fn,K (X, V, p) 



2K+n a -l)/2 n J n J ( nb + 

x (^ /2 pf L^ftH^x) # nt K /2 </) e^-WW)^ (3.16) 

where the H n are Hermite polynomials and the are generalized Laguerre polynomials, n 
indicates the ordered set of vibrational quantum numbers (n B , n a , n^). The eigenvalues are 



e (n,j,K) = ^ + ^ju s + (2n b + if + 1) u h + (n a + u; a + ^ 



J(J + 1)- (^ 2 + ^ 



5. 
(3.17) 



At this order, K is a good quantum number. Thus, we can unambiguously label the 
eigenstates according to n, J, M, and K. The energy levels, 

£~3(£ + e /3), (3.18) 

are labeled by n, J, and K. The operators T^ 4,± \ which enter the analysis at order k 3 , 
will cause the wavefunction to contain contributions from basis functions fij^x 1 w hh K' 
different from the initial choice of K, but each initial choice will correspond to a physically 
distinct energy level. Note that e^' J,K ^ does not depend on a. Higher-order terms will break 
this degeneracy, giving rise to the so-called "/-type doubling" pSfl . 

Perturbation equations for the large-order analysis are obtained from the Schrodinger 
equation by subtituting Eq. (|3.9| ) for the Hamiltonian, and then substituting an expansion 
of the form 

oo 

E = E + k 2 K k e k (3.19) 
for the energy, an expansion of the form 

oo 

$(n,J,K,a) = R: (D-l)/2 g ^(">^)(^ p ) ( 3 . 20 ) 

fc=0 

for the wavefunction, and then collecting terms by powers of k. The zeroth-order solutions 
for the wavefunction are 

^,J,K,a) = ^ (3 21) 

unless K — and a = — , in which case it is zero, according to Eq. Q2.11| ). The e^ ,J ' K "> are 
given by Eq. ( |3.17|) . 

The higher-order results can be obtained recursively. Let us use unnormalized wavefunc- 
tions subject to the condition [p9f] 

(*o|*fc> = h,o- (3.22) 

Then the perturbation equations can be rearranged to obtain a solution for ^/^ in terms of 
the ^j<k- For given choice of n, J, K, and a, let 
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2^ a k Jn',K> 



n'a'.K' 



(3.23) 



The range of the if' summation is from to J, and if if' ^ 0, then a' = + or — . The 
ranges of the n' a depend on the value of k and on the maximum order to which the energy 



expansion is to be calculated [16]. For a calculation up to 62 P , for given p, the wavefunction 
term with the largest number of components is ^/ p , for which n' a < 3p + n a . At order k we 
have 



(if',cr!n / ) 



E 



€k>a k , 



C K i,n> k , =l 

(0) „(*W) , „(+) (K'+l,-v',n") 
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* £1 k',K',n l ,n" u k-k' 



where 



C K ',n' = K - n>s) + (n' a - n a ) u; a + (2n(, + if ' - 2n b - if) u h , 



(3.24) 



(3.25) 



(0) 



'-l) k ^[(k + 3) 7 fc+2 g^ 2 , + ~[j(J + 1) - (if 2 + 13/4)] (k + l)M aK 



+ $k,l 



3A 



~l/2~ 



Qn s ,n'J> 



K,n h ,n'i 



H, 



(±) 



T)*^ s VJ(J + l)-i^±l) 

1/2 



x ^7 fc+1 ti 



M 1 ' . ?(±) _ /^b\ 



1/2 



(3.26a) 



(3.26b) 



in terms of the matrix elements 

p n: n> = 2~ 1/2 (y/n + 15r,+i,„/ - Vri8 n -. 1}n > 

q n>n > = 2~ 1/2 (Vn + 1 <5 n+ i,„' + Vn^n-Ln') 



QK,n,', 



dn.n 1 



3-1/2 
3-1/2 



(2ra + if + 1) £„,„/ - V( n + + + X ) <W 



l.n' 



n(n + if) 5 n -ix 



(2n + K+ 1) <5 nX + J(n + l)(n + if + 1) 5„+i,n' + \/n(n + *0 ^-i,n 



f(+) 
J K,n,' 



(-) 

yK,n,n 

( 7 J 
yK,n,n 



(n + l)(n + 2) <5 n+2 ,n' -y/n(n- 1) <J n -2,n' 
3-1/2 



7i + l)(n + if + 1) 5 n+ l,n' - \M™ + if) <5n-l,n' 

3 _1/2 (v / ^ + if + 15n,n' - V^n-l.n')' 
3 _1/2 (V^+ 15n+l,n' - Vn + if £„,„/), 
3 _1/2 (v / ^ + if + 15 n ,n' + V^^-l/), 
3 _1/2 (v / ^ + 15n+l,n' + + if 5n,n')> 



(3.27a) 
(3.27b) 

(3.27c) 

(3.27d) 

(3.27e) 
(3.27f) 

(3.27g) 
(3.27h) 
(3.27i) 
(3.27j) 
(3.27k) 
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with 7 = {2B/u s ) 1/2 . These expressions are straightforward to derive from the standard 
differential and recursion relations for the Hermite and generalized Laguerre polynomials. It 
is inconsequential that Cx >n = 0, since Eq. ( |3.22|) implies that cl^q' 1 1] = 0. Once Eq. (ggg) 
has been solved up to a given order, the coefficient of the energy expansion corresponding 
to the next higher order can be calculated according to 

c - v V (it® 

e fc — \ u k',K,n,n' a k-k' 
k'=l n' 

4-77-(+) JK+l-a,n') „(-) (K-l-*,n')\ 
~ 1 ~- n k',K,n,n' a k-k' ^ s2 k',K,a,a' a k-k' )• 

(3.28) 

In Appendix B we present a convenient linear algebraic formulation of these recursion rela- 
tions and an analysis of the computational cost. 



IV. RESULTS 
A. C0 2 

An extensive tabulation of rotational spectral lines for CO2 can be obtained from the 
HITRAN data base We will consider here the ground vibrational state and the first 
vibrational (rib = 1). Our potential, Eq. (|3.1|) , depends on the 5 parameters u s , Ub, uj a , R eq , 
and A, with z eq = 0. We will use literature values for the 3 frequencies |JT| and fit R eq and 
A to the empirical spectra. 

We use an iterative fitting procedure, adjusting R cq to fit the splitting between rotational 
levels and A to fit the splitting between vibrational levels. To obtain initial guesses, we 
replace p 2 in the coupling term in Eq. ( |3.1| ) with its diagonal matrix element (2nb+ K + 1) / 'ub, 
which yields 

U(R S ) « ^ s 2 (R s - R nh , K f - \ulRl h ^ (4.1) 



2 



where 



Rn h ,K — Req ~ ~ " 7, — • (4.2) 

w;ojb 

Thus, we have independent constants i?o,o and R^o that can be directly fit to the rib = 
and nb = 1 energies for a given J. These results yield initial approximations for R eq and A, 
and hence and initial calculation for the energy. 

Let R^x and E^ lc (K) be the initial approximations. If we assume that for given J and 
K the energy is proportional to Rq 2 k , then subsequent approximations can be determined 
using 



6(»+l) _ 5(») 



E^ c (K)/E cxpt (K)} 1/2 , (4.3) 
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for K = or 1. For CO2, with J chosen as 84, we found that 3 iterations were sufficient to 
obtain convergence. Our final paramter values are listed in Table |. 

Figures [I] and |2] report our results for the ground and first vibrational states, respectively. 
We show the difference between our calculated energies and energies calculated from a rigid- 
rotor approximation using a rotational constant of 0.39 cm -1 |JT). The empirical results [30 
are indicated by circles. 



B. N 2 and OCS 



We consider here only the ground vibrational states of the linear nonsymmetric molecules 
N 2 and OCS. We set A = and fit the bond distances Ri i3 and B,2,3, where the subscript 
"3" labels the central atom (N or C), and atom 1 is chosen as N for N 2 and O for OCS. 
According to Eqs. (|2.3j), we have 



R C q — ( — - — -) (#1,3 + ^2,3)5 



V m s / 



(4.4) 



-eq 



in 



-Bi, 3 - (— I R23 



mj 



(4.5) 



R 



(i> 



NO 



1.191 A, Rqq = 1.161 A, Rqs = 1.561 A. To fit the rotational spectrum we 

aR$ and obtained the scale factor from 



As initial approximations we used the following bond lengths from the literature [[32 
1.126 A 



p(i) 

%N 



simply scaled these distances according to Rij 
a one-parameter fit to a arbitrarily chosen rotational level. We obtained q;n 2 o 



1.002 and 



«ocs = 1-0256. The resulting values for the R^, which are listed in Table [IJ are within the 
experimental spread of bond length values. For example, our result for Roc agrees with 
the value 1.22 A reported by Hiickel [33]. Figure [I] compares our calculated energies with 
empirical |30| and rigid-rotor |3j| results. 



V. DISCUSSION 

The traditional approach to computing molecular spectra is to numerically diagonalize 
the Hamiltonian matrix in a finite basis set. Our perturbation theory, with an exact recursive 
solution at each order, might appear to have little in common with that method. Note, 
however, that our theory also expresses the wavefunction as a finite linear combination of 
basis functions, according to Eq. ( |3.23|) . Diagonalization uses the variational principle to 
determine the best set of coefficients to multiply the basis functions for a given arbitrary 
finite truncation of the infinite basis set. In contrast, our theory yields the unique asymptotic 



expansion of the wavefunction ||34|| . The finite linear combinations of the basis functions are 



exact solutions for the wavefunction components ^k- Thus, at given order in the perturbation 



theory, the asymptotic principle implies a unique choice of basis set. Chang et al. p5 
incidently, have shown that the basis set chosen by a perturbation theory is superior to a 
minimum-energy basis set for use in numerical diagonalization calculations of vibrational 
spectra. 
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Diagonalization determines the coefficients of the basis expansion that give the "best" 
value for the energy. One can expect that the perturbation theory using the same basis 
set gives a "worse" result |36|. We have found [0] that it is indeed the case, for a model 



vibrational problem, that for a given basis set perturbation theory yields a less accurate 
energy than does diagonalization, and we expect this to be true as well for the problems 
considered here. However, the computational cost of perturbation theory is much lower. 
We show in Appendix B that the cost of calculating the perturbation expansion for a given 
eigenvalue of a triatomic molecule scales as N^ 3 , where N v is the size of the vibrational 
basis. In contrast, the standard diagonalization algorithm scales as iV 3 |]37|| . 

It is of particular interest to consider the effect on cost of the value of J. For variational 
calculations the size of the basis set increases linearly with J ||. Therefore, the cost of 
diagonalization scales as J 3 . For perturbation theory the total computational cost scales only 
linearly with J. Furthermore, the rate of convergence of the perturbation expansion actually 
increases slightly with J P7| . Thus, DPT seems to be an especially appropriate method 
for computing eigenvalues corresponding to high rotational excitation. Pade summation of 
the energy expansion at fifth order yields convergence to 7 decimal digits for J = and 
8 decimal digits for J = 80. We can easily compute these expansions to 15th order on a 
desktop computer workstation, obtaining convergence to more than 12 significant figures. 

The relative cost effectiveness of DPT depends on the magnitude of the accuracy loss, 
compared to diagonalization, for a given basis set. The accuracy of DPT depends on the 
rate of convergence of summation approximants of the perturbation expansion. In contrast 



to dimensional expansions for the electronic Schrodinger equation p8[ , the expansions for 
the systems considered here are rapidly convergent. A comparison of the cost of ob- 
taining a given level of accuracy for vibrational energies of a model 2-coordinate system of 
coupled oscillators showed that, in the end, perturbation theory gave much higher accuracy 
than diagonalization for a given operation count. We have not carried out diagonalization 
calculations for our 3-coordinate systems, so we do not have a detailed comparison of the 
methods in this case. The accuracy vs. cost of the 3-coordinate perturbation expansion is 
comparable to that of the model system, since the rate of convergence for the model system 
is slightly slower but the scaling with N v is slightly better. Therefore, we expect that our 
method is significantly more efficient than straightforward diagonalization in this case as 
well. 

The three molecules treated here are tightly bound and therefore can be accurately 
described with harmonic, or nearly harmonic, potentials. The rate of convergence for po- 
tential surfaces with larger anharmonicity will undoubtedly be slower. In general one can 
expect, from fundamental theoretical considerations PS,3S|, that dimensional expansions 



and coupling-constant expansions for anharmonic oscillators have a radius of convergence 
of zero. Coupling-constant expansions calculated by Cizek et al. || for model vibrational 
problems with somewhat larger anharmonicity did indeed exhibit only semiconvergent be- 
havior, with divergence setting in at higher orders, although the expansions were rapidly 
summable with Pade approximants. With truly large anharmonicity as in van der Waals 
molecules the divergence will undoubtedly be more of a problem. However, there are tech- 
niques that could be used to improve the rate of convergence in such cases. For example, 
one could redefine the dimensional continuation so that some dependence on J is present 
at zeroth order. This would make it possible to include centrifugal distortion in the basis 
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functions. Another useful technique is to modify the potential energy operator so as to 
include an arbitrary renormalization parameter, the value of which is chosen to optimize the 
convergence HOI. In the context of DPT, this consists of using a dimensional continuation of 



the Hamiltonian that has explicit dimension dependence in the potential, giving the correct 
potential at D = 3 but increasing the stability of the system at large D [^5|] . 

In recent years, methods have been developed |3|,|5|, fni , f^ , ft3f| that lower the cost of di- 
agonalization for molecular spectra, especially in the case of floppy molecules, for which a 
harmonic-oscillator basis set converges very slowly. In view of the similarities between the 
form of the wavefunction in perturbation theory and in diagonalization, it may be possible to 
make analogous improvements in the perturbation theory. In any case, perturbation theory 
has the advantage of simplicity and directness, in that a straightforward recursion relation 
yields the exact values of the coefficients of the unique asymptotic expansion. Furthermore, 
no modification of the theory is needed in order to study resonances, in which case quadratic 
Pade summation of the perturbation expansion yields complex energies |7|Ji^]. 

Note also that, in contrast to matrix diagonalization, the perturbation theory provides 
each eigenvalue with an unambiguous label, in terms of the quantum numbers (n, J, K, a). 
This labeling is analogous to that provided by the vibrational SCF method ||45|| ; however, 
that method introduces a separability approximation while the perturbation theory, at least 
in principle, will converge to the exact result. In practice, the theory presented here can 
break down if the zeroth-order, harmonic, eigenstate is degenerate or nearly degenerate 
with another harmonic eigenstate [}46| |. Such cases can be treated by including additional 



dimension dependence in the potential energy so as to break the degeneracy in the large- .D 
limit, although this procedure results in a somewhat slower rate of convergence. 
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APPENDIX A: KINETIC ENERGY OPERATOR IN BODY-FIXED 

COORDINATES 

In order to construct the basis functions 4>j,m,k{Rs, P, z), it is convenient to introduce 
Wigner functions P?| |, 



4m'(M,7) = (L,M|e^V^V^|L,M'>, (Al) 



d L MM m = e ~ tMa v m, m <{^ & ^* e-lM ' 7 ' ( A2 ) 

where (a, /3, 7) are Euler angles and the Li are components of the orbital angular momentum 
operator L along perpendicular axes Xj. The T> functions are eigenfunctions of L and Lj, 
and of the component Ly along a body-fixed axis x^, with eigenvalues L(L + 1), M, and 
M', respectively. Let (Pj,m,k(R s , z, p) be a set of orientation-independent basis functions, 
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$(Rs, Ra) 



9 7 4- 1 \ 1 / 2 

" " T, V M,K^Th,M,A^^P)- (A3) 
if' 



The kinetic energy operator in this basis is 

2J+ 1 



Tj,k,k> = j v J M>K (n) r rcl v J MjK ,(ny dn, (A4) 



which is related to T re \ according to 
'2J+ 1^ 1/2 



8tt 2 



v J MjK (n)T iel ^(R^,-R a ) dn 

= E T J,K,K> 4>J,M,K'{Rs, Z, P)- (A5) 
K' 



The wavefunction can also be expanded in terms of a set of basis functions 
0j ) M,L s ,L a (-R s , -R a ), with quantum numbers L s and L a for the angular momenta corresponding 
to the vectors R s and R a . Let fl s and f2 a be the Euler angles that describe the orientation 
of R s and R a , respectively, relative to the space- fixed axis. Then 



$(R S , Ra) = — E V( 2L s + X )( 2L - + X ) 0JM,L s ,L,(Rs, 
47F J,M,L B ,L a 



E fiwtW zfc oW (A6) 



x 

M s ,M a 

where the C^ sMs . La Ma are Clebsch-Gordan coefficients. 

The angles that describe the orientation of R a relative to R s are fi a — fl s = (a, /?, 0) as 
illustrated in Fig. [| The polar angle j3 is the angle between R a and R s . Using the addition 
theorem for T> functions and the Clebsch-Gordan series for the product of D functions, we 
find that 

$(R s ,R a ) = E V2J+lV J M}K (a s ,P s ,a)* 

47F J,M,K 



x E (-l) La+X V2L a + 1 C L jf K .^ K d L K %(5) 4>j,m,lMRs, J?0 • (A7) 
Therefore, we can express the former basis functions in terms of the latter according to 



K (R S , z,p)= E {-l) L * +K \JL a + 1/2 ("rV.L .K d L K *M Oj.m.^JH. R & ). (A8) 



On the left-hand side of Eq. (|A8j) we have replaced the internal coordinates _R a and f3 with 
the cylindrical coordinates defined in Eqs. ( [2.10| ). 

Equation (|A8[ ) is useful because the representation of T re i with respect to the 0j,M,L s ,L a 
is particularly simple. Suppose we write 

rel = ~2RldR~ s dR s + 2B* 8 " 2J%M~ a a dR a + 2Rj a ( } 
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and substitute this into Eq. (|A4j) . Then if we substitute Eq. ( |A8|) into Eq. ( |A5|) , we can 
replace the operator L 2 with L S (L S + 1). It follows that 

/2J+1\ 1/2 



V 8tt 2 



E 



P^(0)L s 2 $(R s ,R a )rff2 
(2L s + l)(2L a + l)l ! 2 



2(2j + 1; 



L S (L S + 1) Ci>" K d£M 0J,M,L s ,L a (R S , Ra). (A10) 



This result is in terms of the 0j,M,L s ,L a > but these functions can be expressed in terms of the 
0j,m,k by inverting Eq. (TA"8~D, 



>J,M,Ls,L, 



(Rg, R a 



L a + l/2 ^(-l) La+K ' C is '° 



J,-K';L & ,K' 



x / <t>J,M,K'{Rsi -RaCos/5, i? a sin P)d£, (j3)fanPdp, (All) 



which follows from the orthogonality condition for the X> functions |48|| . 

We now evaluate the sum over L s that results from substituting Eq. ( |A11| ) into Eq. ( |A10| ). 
This substitution gives a product of the Clebsch-Gordan coefficients can be expressed as an 
integral over T> functions, 



-1) J 



Note that 



2L S + l ^ 1 / 2 JK Ls0 

2 J + l J L 'L s ,0;L„K^ J-K';L a ,K 



( 1) ^ 'J,-K;L*,K^ 'J,-K';L a ,K' 

J v^(ny v J _ K _ K ,(n') p£>(n') do'. (Ai2) 



L S (L S + l)Vfc = L 2 s Vfa = (Lq - L +1 L_! - L-rL+i P, 



0,0 ' 



(A13) 



where Lo and L±i are the components of L s in terms of spherical coordinates. It follows 
from Eqs. (|A12|) and (|A13|) , and the orthogonality condition, that 



E(-i 



, /2L S + 1\ 1/2 



2 j i j L s (L s + 1) C L ' B)0;L&)K Cfl K ,. L& K , — ^2 L S (L S + 1) Cj*l K . L& K Cj a l K ,. L& K , 

[J(J + l)-2K 2 + L a (L a + l)]<^ 
-{[J( J + 1) - K{K - l)][L a (L a + 1) - - 1)]} 1/2 
-{[ J( J + 1) - K(K + 1)] [L a (L a + 1) - K(K + 1)]} 1/2 fo, )JC+1 . 



(A14) 



Thus, the term L 2 /(2R 2 ) in T rel couples basis functions with neighboring values of K. 

Next we need to remove the quantum number L a from the expression for Tj,k,K'- We 
deal with the operator LJ (2i? 2 ) by simply replacing it with its differential form, 



d 2 



da 2 



2 <9 2 2 (I d d Id 2 



dz 2 



p dp dp p 2 da 2 



2 4 + {» 



d_ 

dp 



1+1 



D L K fi (a,f3y. 
(A15) 
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To remove the L a dependence that comes from the operator L^/(2R^), via Eq. ( |A14j ), we 
use the fact that 



>L a (L a +l)-K(K±l) DfrJa, /3)* = ±2^ L Tl D^ ±1 Ja, $)* 



±e 



P 



d_ 

dz 







z d 



dp pda 



e i(K±i )a d L a± M ] (A16) 



The integration over (3, from Eq. ( |A11| ), is evaluated using the completeness property of the 
d functions: 



£ (L + 1/2) dZxjoip) 4k,o(/3') = $(cosf3 - cos/3')- 



(A17) 



L=K 



We use the fact that i(d/da)V^ = KV^ to evaluate the a derivatives. 

The two remaining operators in Eq. ( |A~9| ) are simple to evaluate in the new basis. The 
derivatives with respect to R s are the same in either basis and the derivatives with respect 
to i? a are independent of the orientation of the body-fixed coordinate system. We can make 
the substitution 



1 R2 



2RldR a 



dR f 



+ 



1 d 2 



1 d d 



1 d 2 



2RI 2dz 2 2pdp P dp 2p 2 da 2 



(A18) 



Eq. ( 2.12) ) follows after the transformation to the symmetrized basis functions of Eq. ( ^.11 ) 



APPENDIX B: LINEAR ALGEBRAIC FORMULATION OF THE RECURSION 

RELATIONS 



A linear algebraic method for solving perturbation equations was recently presented 
by Dunn et al. [[L(J. This approach is completely equivalent to Eqs. ( |3.24| ) and ( |3.28| ). 
However, it expresses these equations in a transparent form that simplifies the analysis of 
computational cost and is particularly convenient for efficient implementation on a computer 



1. The matrix method 



This method uses the fact that each can be exactly expressed as a finite linear 
combination of the eigenfunctions of 7io, according to Eq. ( |3.23j ). Let aj^ ,<T ^ be the rank-3 
tensor of expansion coefficients a\^' a ' n \ with the tensor index (ii,i2,is) corresponding to 
n' = (ii — 1,22 — 1)^3 — 1)- Then operators can be expressed in terms of matrices. In 
particular, Eqs.( p.24| ) and ( |3.28| ) can be written as 



'K> 



k'=l 



Cfc'i ® i <S> i — H 



(o) 

k'K> 



,(*V') 
l k-k' 



tt(+) (K'+l,-*') 



l k'K' 



Hi./ 



(JC'-l, 
k-k' 



V) 



(Bl) 
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V fH (0) a (K,CT) + H {+) a iK+1 '- a) + H H a^ 1 '^ fB2l 

^fe'^^-fe' + n k',K A k-k' + "-k^K^-k-k' ) ' l OZ J 



fc'=l 



where 



H 



(o) 

k,K 



-l) fc ^ s7 2 ^ (k + 3) 7 fc+2 q fe+2 g) i ® i + ^ J(J + 1) - (K 2 + 13/4) (jfe + l) 7 fe q* ® i ® i 



+4/c 7 



fe-i 



i! * q fc ~ x <g> p 2 ® b x - ( — J q^" 1 ® z 2 ® c x + q fc ~ x ® d ® 



H 



(±) 



3A • i~ 

+ *m i/2- q ® i ® b^, 

-l) fc ^ sv /J(J+l)-^±l) 



x^7 



fc+i 



Wb/ \u; a / 



1/2 



o (±) 



(B3) 



(B4) 



Ck> = E [« - n a )oj s + « - ™ a R + (2n b + AT' - 2n b - X)a) b ] (8) i«) ® i(n b ) 



(B5) 

with 7 = (2B /ujs) 1 ^ 2 , in terms of the matrices [i]^- = Sij and = ^n+i^j, and the 

matrices p, q, hx, ck, d, ex, , Sk \ an d z defined by Eqs. (|3.27|) . We use the direct 
product notation x s ® x a ® x b a to represent the following procedure: multiply the columns 
of a, which correspond to the quantum number n b , by the matrix x b ; multiply the west-east 
rows, corresponding to n a , by x a ; and then multiply the north-south rows, corresponding to 
n s , by x s . The dot-product notation in Eq. ( p2|) represents a scalar product defined by 



The dependence of the operators Hjj?^ and Jl^K is quite simple. It is found only in the 
operations by the coordinate x, on account of the way in which we defined the dimensional 
continuation of the kinetic energy operator. We can take advantage of this to streamline the 
form of the recursion relations. Let us define the operators 



,(o) 



1 

-L 

2 



,,3„3 



1+ 3 7 



J(J + 1) 



K 2 H 

4 



q® i ® i 



+ 



i ® p 



•K 



h (±) 



i< ---^J(J+l)-K(K±l) 1 



i (g> z <g> c K + i ® d (g) 

1/2 / ,7), \ V2 



and the tensors 



1 „ f o .4 
--^ s7 ^|3 7 d q®i® i + - 7 



(+)(K+l,-a) 



^a 

— J 1 » P 



, (-) (*-l-«r) 



fM - [ - 



w b 



J(J+1)- K 2 + 
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q <g) l ® l > a] 



(B7) 
, (B8) 

(B9) 
(BIO) 
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Then the recursion relations can be written as 



I j=0 J 



(Bll) 



and 



AK,a) 



{ ® i ® bi, a£f - £ [(* - i) Ag' CT) + Bg CT) ] ) , (B12) 

LcUs CUb j=0 J 



in terms of the tensors 



A<^> = q*-'- 1 ® i ® i af y \ 
Bg^ = q^- 1 ®i®ibf f ^ I 



(B13a) 
(B13b) 



where q = — 7q. Equations ( |B 1 1|) and ( B12|) are convenient because the Ajg V) and Bjg V) 
can be calculated using recursion relations that are particularly simple: 



A^ = q<8.i®iA^, 



B^ = q®i®iB£V). 



(B14a) 
(B14b) 



2. Computational cost 

A close examination of Eq. (|B1|) reveals that the columns of aj^ ,<J ^ will have one more 
nonzero element than will the columns of a^_^\ the west-east rows will have two more 
nonzero elements, and the north-south rows will have three more nonzero elements. There- 
fore, the total number of nonzero elements in a^ scales with order k as k 3 . It follows 
that the number of multiplications needed at order k in the recursion to compute e\ \ 

or to compute A^'' 7 ^ and B^ ,<J ^ using Eqs. ( P14j) , scales as A; 3 . 

Consider the cost of evaluating Eq. (|B11| ) for all k such that k < p, where p is the order 
at which the energy expansion will be truncated. The most expensive part of the calculation 
is the evaluation of the 3 terms that are summed over j. For each value of j there will be 
a single matrix operation, scaling at worst as k 3 , and the sum over j is repeated for each 
value of k. Summing the cost over all the values of j and k gives a cost scaling of p 5 for 
each value of (K',a'). The range of K' is < K' < J. Therefore, the total cost scales as 
Jp 5 . According to Eq. ( |3.23| ), the number of basis functions, N v , is equal to the number 
of nonzero elements of the corresponding tensors. Therefore, N v scales as p 3 , and the total 
cost scales as JN^J 3 . 
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TABLES 



TABLE I. Frequencies, bond distances, and anharmonicity constant. 



molecule 


uj s (cm x ) 


(cm x ) 


to h (cm x ) 


#1,3 (A) 


#2,3 (A) 


10" 6 x A (cm" 5 / 2 ) 


C0 2 


1351. a 


2396.4 a 


672.2 a 


1.163647 


1.163647 


1.79959 


N 2 


1285.0 b 


2223. 5 b 


588.8 b 


1.128249 


1.193379 





OCS 


859.0 b 


2079.0 b 


527.0 b 


1.190693 


1.600922 






Ref. 0. b Ref. |2|. 
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FIGURES 



FIG. 1. Deviation from rigid-rotor approximation of the energy levels of the ground vibra- 
tional state, vs. the angular-momentum quantum number J. The ordinate corresponds to 
E(J) - E(0) - B Iot J{J + 1), with £ rot (C0 2 ) = 0.3906 fl[, B ro t(N 2 0) = 0.419 @, and 
-Brot(OCS) = 0.20286 |3^|. The solid curves show our calculated results while the circles show 
empirical results [30]. 



FIG. 2. Deviation from rigid-rotor approximation of the energy levels of the first excited vi- 
brational state, vs. the angular-momentum quantum number J. The ordinate corresponds to 
E(J) - E(0) - B TOt J(J + 1), with B rot (C0 2 ) = 0.3906 @. The solid curves show our calculated 
results while the circles show empirical results p0[. 



FIG. 3. Euler angles (a s , f3 s ) for the orientation of the body-fixed axis R s relative to space-fixed 
Cartesian axes, and {a, (5) for the orientation of the vector R a relative to R s . 
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